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Abstract 


The  paper  describes  a multiple  grid 
method  for  transonic  flow  calculations. 

The  proposed  scheme  incorporates  a general- 
ized alternating  direction  method  as  the 
smoothing  algorithm.  Numerical  experiments 
with  this  multigrid  alternating  direction 
(MAD)  method  indicate  that  it  is  both  fast 
and  reliable. 


1.  Introduction 


The  multiple  grid  method  was  first 
proposed  by  Federenko  (1],  who  realized 
that  it  should  be  possible  to  accelerate 
an  iterative  scheme  for  solving  difference 
equations  by  calculating  corrections  for 
the  fine  grid  equations  on  a sequence  of 
successively  coarser  grids.  This  idea  was 
subsequently  analyzed  by  Bakhavlov  [2], 
and  then  extended  and  applied  to  a variety 
of  problems  by  Brandt  [3].  It  has  recently 
been  proved  under  rather  general  assump- 
tions by  Nicolaides  [4],  and  Hackbusch  [5], 
that  the  number  of  operations  required  to 
solve  the  equations  arising  from  a finite 
element  or  a finite  difference  approxima- 
tion to  an  elliptic  problem  by  a multiple 
grid  method  is  directly  proportional  to 
the  number  of  unknowns. 


on  fine  grids  before  passing  to  coarser 
grids , and  South  and  Brandt  were  able  to 
obtain  convergence  for  a wider  range  of 
cases  by  using  multiple  line  relaxation 
sweeps  in  different  directions  [7]. 

In  the  present  work  the  difficulties 
have  been  attacked  by  combining  the  multi- 
ple grid  method  with  a generalized  alternat- 
ing direction  method  suitable  for  transonic 
flows  as  the  smoothing  algorithm.  Numerical 
experiments  indicate  that  this  multigrid 
alternating  direction  (MAD)  method  converges 
rapidly  and  reliably  for  a range  of  cases 
typical  of  the  cruising  regime,  up  to  the 
onset  of  drag  rise.  It  appears  also  that 
the  method  can  readily  be  generalized  to 
treat  three  dimensional  flows. 

2.  The  Difference  Approximation 

The  case  which  will  be  considered  is 
that  of  two  dimensional  transonic  flow 
past  an  airfoil,  using  the  potential  flow 
approximation,  which  has  been  found  to  give 
useful  predictions  in  practice  for  flow 
containing  shock  waves  of  moderate  strength 
[8].  -The  potential  flow  equation  will  be 
treated  in  the  conservation  form 

(pu)  + (pv)  - 0 (1) 


There  is  less  experience  of  the  use  of 
multiple  grid  methods  for  nonelliptic  prob- 
lems. The  first  demonstration  of  the  use 
of  a multiple  grid  method  for  a transonic 
flow  problem  was  by  South  and  Brandt  16) , 
who  solved  the  transonic  small  disturbance 
equation  for  a nonlifting  flow  and 
observed  a high  rate  of  convergence. 
Difficulties  were  experienced,  however, 
both  by  South  and  Brandt  and  by  the 
present  autnor,  in  the  treatment  of  lifting 
flows  and  in  calculations  on  nonuniform 
and  curvilinear  meshes.  There  was  a 
tendency  to  produce  an  oscillating  sonic 
line,  and  for  the  calculations  to  enter  a 
variety  of  limit  cycles  between  several 
grids.  These  difficulties  appeared  to  be 
due  to  insufficient  smoothing  of  the  errors 
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where  x and  y are  Cartesian  coordinates, 
p is  the  density,  and  the  velocity  compon- 
ents u and  v are  the  gradient  of  the 
potential  4, 


u - *x  * v - 4 . (2> 

If  q is  the  speed  /u2+v^,  the  local 
speed  of  sound  a is  determined  by 
Bernoulli's  equation 

-2  - «}  - ^ q*  ‘3) 

where  aQ  is  the  stagnation  speed  of  sound. 
The  density  follows  from  the  relation 


,T-1 


Mia2 


(4) 


where  M_  is  the  Mach  number  q/a  of  the 
uniform  flow  at  infinity  and  y is  the  ratio 
of  specific  heats.  At  tl.e  profile  the 
potential  satisfies  the  Neumann  boundary 
condition 


>4 
5 n 


0, 


(5) 
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where  n is  the  normal  direction,  and  also 
ths  Kutta  condition  that  the  tangential 
velocity  is  bounded  at  the  trailing  edge. 

At  infinity,  the  potential  approaches  the 
potential  of  a vortex  in  compressible  flow. 

In  the  numerical  experiments  reported 
here  these  equations  are  solved  in  a 
coordinate  system  generated  by  a conformal 
mapping  of  the  domain  onto  the  interior  of 
a circle.  Using  polar  coordinates  r and  6 
the  potential  flow  equation  becomes 

fe  <P*9>  ♦ r |?  " 0 (6) 


The  velocity  components  in  the  6 and  r 
directions  are  , 

r$fl  r ♦ 

u - -J-  , v - — (7) 


H 


where  H is  the  modulus  of  the  transforma- 
tion onto  the  exterior  of  the  circle. 


The  difference  approximation  is  similar 
to  schemes  which  have  been  previously  used 
[P).  It  is  derived  by  augmenting  a central 
difference  scheme  with  an  artificial 
viscosity  which  introduces  an  upwind  bias 
throughout  the  supersonic  zone.  Using  sub- 
scripts i.j  to  denote  values  at  mesh  points, 
and  i+1/2,  j+1/2  to  denote  values  at  the 
midpoints  of  the  segments  connecting  mesh- 
points,  the  approximation  to  equation  (6) 
has  the  form 


{pi+l/2,j(*i+l,j-*i,j) 

‘ pi-l/2,j(*i,j‘<,i-l,j)} 

+ ^ {rj+l/2pi,j+l/2(*i,j+r*i,j> 

" rj-l/2Pi,j-l/2(*i,jH’i,j-l)}+  Tij“  0 (8) 

where  a 9 and  &r  are  the  mesh  widths,  and 
T..  is  the  artificial  viscosity.  The  flow 
i.Vthe  supersonic  zone  is  assumed  to  be 
roughly  in  the  6 direction  and  an  artifi- 
cial viscosity  which  gives  a bias  only  in 
the  6 direction  has  been  used  in  the 
experiments  so  far.  Let  u be  a switching 
function 

u - max  (0,  (l-M2/M2) } (9) 

which  vanishes  when  the  local  Mach  number 
M is  below  a cutoff  Mach  number  Mc  , and 
let  S-  , be  i difference  approximation  to 
1 » J 
2 2 

uc (u  /a  )«90  at  the  point  i.j.  Then  the 


artificial 

viscosity  is  of 

the  form 

Ti.j 

where 

“ Pi*l/J,j  ’ Pi 

“1/2 , j <10> 

, -1 

Si.j  "tSi-l,j 

if  ui+l/2,j>0 

Pi+l/2,j  | 

Sltl,J_cSi+2.j 

if  Ui*l/2,J<0 
(11) 

r | « m _r>» ram«f  COr.trcliir.j  2CCU? 

acy.  If  c ■ 0 the  added  terms  are  of  order 
16,  yielding  a first  order  accurate  scheme, 
and  if  e • 1 the  added  terms  are  of  order 
19^.  yielding  a second  order  errnrere 
scheme.  Inc  cutoff  value  M-  » 0.9  nas  been 
found  to  result  in  reliablecconvergence  of 
the  multigrid  alternating  direction  scheme. 

3.  Review  of  the  Multiple  Grid  Method 


Consider  the  solution  of  the  equation 

Lhu  - f (12) 

by  a relaxation  method,  where  Lh  approxi- 
mates a linear  differential  operator  L on  a 
grid  with  a spacing  proportional  to  the 
parameter  h.  Let  U be  an  approximation  to 
the  solution,  and  let  V be  the  correction 
to  U such  that  U ♦ V satisfies  (12).  Then 
the  basis  of  the  multiple  grid  method  is  to 
replace  (12)  by 

L2hV  + I^L1^  - f (13) 

n 

2h 

where  L is  the  same  approximation  to  L on 
a grid  in  which  the  spacing  has  been 
doubled,  and  If-*1  is  an  operator  which 
transfers  to  each  grid  point  of  the  coarse 
grid  the  residual  LhU  - f of  the  coincident 
point  of  the  fine  mesh,  or  alternatively 
a weighted  average  of  the  residuals  at 
neighboring  points.  After  solution  of  (13), 
the  approximation  on  the  fine  grid  is 
updated  by  interpolating  the  correction 
calculated  on  the  coarse  grid  to  the  fine 
grid,  so  that  U is  replaced  by 

Un*w  - C ♦ I^V  (14) 

where  I»h  is  an  interpolation  operator. 
Equation  (13)  can  in  turn  be  solved  by 
introducing  an  approximation  on  a yet 
coarser  grid,  so  that  a multiple  sequence 
of  grids  may  be  used,  leading  to  a rapid 
solution  procedure  for  two  reasons.  First, 
the  number  of  operations  required  for  a 
relaxation  sweep  on  one  of  the  coarse  grids 
is  much  smaller  than  the  number  required  on 
the  fine  grid.  Second,  the  rate  of  conver- 
gence is  faster  on  a coarse  grid,  reflect- 
ing the  fact  that  corrections  can  be 
propagated  from  one  end  of  the  grid  to  the 
other  in  a small  number  of  steps. 

To  extend  this  idea  to  nonlinear  equa- 
tions, equation  (13)  may  be  reorganized  by 
adding  arid  subtracting  the  current  solution 
U to  give 

L2h(U  «■  V)  ♦ I^lN)  - L2hU  • f 
or  ,hh  , 

L2hff  - i (15) 

where  5 is  the  improved  estimate  Of  the 
solution  to  be  determined  on  the  coarse 
grid,  and  1 is  an  appropriately  modified 
right-hand  side, 

1 - f ♦ L2hU  - I^lNj  (16) 

The  updating  formula  (14)  now  becosws 
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Unew  - 0 + 1^(0  " u)  (17) 
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perturoation  operator  to  represent  the 
correction  equation  (13) • 


4.  Smoothing  Algorithms 

The  success  of  the  multiple  grid  method 
generally  depends  on  the  use  of  a relaxa- 
tion algorithm  which  rapidly  reduces 
the  high  frequency  components  of  error  on 
any  given  grid,  because  on  a coarser  grid 
these  components  cannot  be  distinguished 
from  low  frequency  components.  This  alias- 
ing process  will  cause  improper  correc- 
tions to  be  computed  on  coarse  grids,  and 
can  prevent  convergence. 


It  turns  out  that  point  and  line  relax- 
ation schemes  do  not  necessarily  provide  the 
required  smoothing  of  all  high  frequency 
components  of  error  on  a nonuniform  or 
curvilinear  mesh.  To  illustrate  this 
consider  the  model  problem 


a*xx  + b\y  ' 0 (l8) 

with  positive  coefficients  a > 0,  b > 0. 

Let  ££  and  Sx  denote  second  difference 
operators  in* the  x and  y directions,  and 
suppose  that  the  difference  approximation 
has  the  form 

M 2 (Ai*  + B - 0 (19) 

x y 

where  if  Ax  and  Ay  are  the  mesh  widths , 

A “ -■* k , B * — — y (20) 

Ax‘  Ay 

Assuming  periodic  boundary  conditions 
suppose  that  the  solution  after  n itera- 
tions has  the  form 


4 - Gn  eipx  eiqy 


(21) 


where  G is  the  amplification  factor,  and 
let 


P Ax  - 5 


q Ay  ■ n 


Then  a Gauss  Seidel  scheme  yields 
Aei5  + Bein 

G A(2-e‘i5)  ♦ B(2-e'in) 


Suppose  that  Ax  , Ay  are  such  that  A >>  B 
and  consider  the  case  of  a high  frequency 
in  the  y direction  and  a low  frequency  in 
the  x direction,  which  may  be  represented 
by  ( ■ 0 , n • it.  Then 


G 


A - B 
A + 3B 


-v  1 . 


A similar  analysis  of  a line  relaxation 
scheme  shows  that  if  A >>  B effective  damp- 
ing of  high  frequency  error  components 
requires  the  use  of  a scheme  solving  along 
the  lines  in  the  x direction  [6]. 


On  a nonuniform  grid  A and  B may  vary 


widely  so  that  in  some  regions  A >»  B and 
in  others  B >>  A.  This  has  led  South  and 
Brandt  to  use  multiple  line  relaxation 
sweeps  in  different  directions  £71.  An 
alternative  approach  proposed  by  Arlinger 
[9]  is  to  use  auxiliary  grids  constructed 
by  increasing  the  mesh  width  in  one  direc- 
tion only  while  the  line  relaxation 
scheme  is  applied  to  the  lines  in  the  other 
direction.  This  method  has  been  tested  by 
the  present  writer  and  found  to  give  a use- 
ful acceleration  of  the  rate  of  convergence 
of  transonic  flow  calculations.  Its  poten- 
tial efficiency,  however,  is  less  than  that 
of  a full  multigrid  scheme  in  which  the 
mesh  width  is  increased  in  both  directions, 
because  the  coarse  grids  contain  more  mesh 
points  when  the  mesh  width  is  only 
increased  in  one  direction. 

It  is  proposed  here  to  use  an  alternat- 
ing direction  method  as  the  smoothing 
algorithm.  Consider  the  model  equation 
(18)  and  suppose  that  the  correction  <4  is 
calculated  by  the  equation 

(a  - Afi*)  (a  * Bfi*)  ««  “ woLC  (22) 
x y 

where  a is  a parameter  to  be  chosen,  u is 
an  overrelaxation  factor,  and  the  residual 
LA  is  calculated  using  the  result  of  the 
previous  iteration.  (The  usual  Peaceman 
Rachford  scheme  110]  is  obtained  by  setting 
u - 2.)  Then  on  inserting  a trial  solution 
of  the  form  (21)  it  is  found  that  all  high 
frequency  components  are  rapidly  damped  if 

0^4  m(n(A,8)  , w ' 1.5 

(assuming  that  A > 0,  B > 0). 

5.  Generalized  Alternating  Direction  Scheme 

In  the  case  of  transonic  flow  we  have 
to  allow  for  a change  from  elliptic  to 
hyperbolic  type  as  the  flow  becomes  locally 
supersonic.  In  the  model  problem  (18)  this 
corresponds  to  one  of  the  coefficients, 
a say,  becoming  negative.  The  classical 
alternating  direction  scheme  (22)  then  has 
the  disadvantage  that  if  one  regards  the 
iterations  as  representing  time  steps  At  in 
an  artificial  time  direction  t 111],  it 
simulates  the  time  dependent  equation 

O At  - a + b 4yy 

When  a < 0 and  Cauchy  data  is  given  at  x«0, 
corresponding  to  supersonic  inflow,  this 
leads  to  an  ill  posed  problem  which  admits 
oscillatory  solutions  which  are  undamped 
in  time  and  grow  in  the  x direction. 

The  following  generalised  alternating 
direction  scheme  is  therefore  proposed. 

Let  the  scalar  parameter  o in  equation  (22) 
be  replaced  by  a difference  operator 

S i a.  ♦ a,  fi”  ♦ a.fi”  123) 

0 1 x 2 y 

where  and  !"  denote  one  sided  difference 
operators  in  "the  x and  y directions. This 
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yields  the  scheme 

(S-AJ*)  (S-B«5>  «♦  - u S t ♦ (24) 

x y 

in-  which  the  residual  -L«  i»  uiricieuCe  «S  by 
the  operator  S.  The  corresponding  time 
dependent  equation  is  now  a hyperbolic 
equation  of  the  form 

9.3  + 3,3  . + g-3  . m ad  + bd 

0 t lyxt  2vyt  wxx  *yy 

where  the  coefficients  B.tSj/S,  depend  on 
the  parameters  Og.a^.aj.0  1 i 

The  artificial  viscosity  introduced  in 
Section  2 is  equivalent  to  a switch  to 
upwind  differencing  in  the  x direction  as 
was  first  proposed  by  Murman  and  Cole  [12], 
and  in  the  implementation  of  the  scheme 
the  difference  operator  in  the  first 
factor  is  correspondingly  replaced  by  an 
upwind  difference  operator  when  A < 0,  and 
the  direction  of  the  one  sided  operator  5~ 
in  S is  chosen  to  be  upwind. 

The  generalized  scheme  (24)  can  be 
related  to  previously  proposed  alternating 
direction  methods  for  transonic  flows  [13, 
14]  by  appropriate  specialization  of  the 
parameters.  For  example,  the  AF2  scheme 
proposed  by  Ballhaus,  Jameson  and  Albert 
[13],  corresponds  to  the  choice  a.  » 0, 
a,  ■ 0,  followed  by  an  integration  in  the 
x^direction  which  eliminates  the  differenc- 
ing of  the  residual. 

6.  Multiple  Grid  Strategy 

Theoretical  estimates  of  the  rate  of 
convergence  attainable  by  the  use  of  multi- 
ple grids  have  been  obtained  for  recursive 
strategies  [4,5].  South  and  Brandt  [6] 
used  an  adaptive  strategy,  with  transitions 
to  a coarser  grid  if  the  rate  of  conver- 
gence becomes  too  low  on  a particular  grid, 
or  to  a finer  grid  if  the  average  residual 
has  been  sufficiently  reduced. 

In  the  present  work  a simple  fixed 
strategy  has  been  found  to  be  effective. 
Each  cycle  begins  on  the  fine  grid.  The 
alternating  direction  iteration  is 
performed  once  on  each  grid  until  the 
coarsest  grid  is  reached.  Then  it  is 
performed  once  on  each  grid  going  back  up 
to  the  second  finest  grid,  and  the  cycle 
terminates  with  the  interpolation  of  the 
correction  from  the  second  finest  grid  to 
the  fine  grid.  It  is  convenient  to  measure 
the  work  in  units  representing  the  work 
required  to  perform  one  iteration  of  the 
alternating  direction  scheme  on  the  fine 
grid.  Since  each  grid  has  1/4  as  many 
cells  as  the  next  finer  grid,  the  work 
required  to  perform  each  cycle  is 

1 + 2(?  * I?  * ST  **•]  i unit*' 

plus  the  overhead  of  computing  and  trans- 
mitting residuals  from  one  grid  to  the 
next,  and  interpolating  the  corrections. 

It  is  the  usual  practice  to  accelerate 


the  alternating  direction  scheme  (22)  by 
using  a sequence  of  values  of  the  parameter 
a designed  to  give  rapid  damping  of  the 
error  components  in  a series  of  frequency 
w&nda • The  multrgrid  alternating  direction 
method  economizes  the  work  required  by 
passing  to  the  coarse  grids  to  treat  the 
error  components  in  tha  low  frequency 
bands.  If  a sequence  of  6 parameters  were 
used  to  treat  6 frequency  bands,  for 
example,  the  work  required  to  complete  one 
cycle  through  the  parameters  would  be 
6 units,  whereas  the  work  required  to 
perform  the  alternating  direction  itera- 
tions of  a multigrid  cycle  with  6 grids 
would  be  less  than  1 2/3  units  with  the 
present  strategy. 

7.  Results 

The  efficiency  of  the  multigrid  alter- 
nating direction  method  has  been  confirmed 
by  numerical  experiments.  Some  typical 
results  are  presented  here.  All  the 
examples  were  calculated  on  a circular 
domain  generated  by  conformal  mapping  of 
the  profile  to  a unit  circle,  with  192 
cells  in  the  9 and  32  cells  in  the  r direc- 
tion on  the  fine  grid.  Up  to  6 grids  were 
used,  giving  a coarse  grid  with  as  few  as  6 
cells  in  the  9 direction  and  1 cell  in  the 
r direction. 

Figure  1 shows  the  observed  convergence 
rate  using  different  numbers  of  grids  for  a 
case  with  a shock  of  moderate  strength,  the 
flow  past  an  NACA  64A410  airfoil  at  Mach 
0.720  and  an  angle  of  attack  of  0*.  The 
pressure  distribution  of  the  fully 
converged  result  is  shown  in  Figure  la. 
Figures  lb-lg  show  the  logarithm  of  the 
average  absolute  value  of  the  residual 
plotted  against  the  work,  measured  by  the 
equivalent  number  of  iterations  on  the  fine 
grid.  The  convergence  rate,  measured  as 
the  mean  reduction  in  the  average  residual 
oer  unit  of  work,  is  also  indicated  under 
each  graph.  It  can  be  seen  that  the  rate 
of  convergence  was  improved  from  0.9840  to 
0.6677  as  the  number  of  grids  was  increased 
from  1 to  5,  while  no  further  improvement 
was  realized  with  6 grids.  (In  subsonic 
flow  it  pays  to  use  6 grids.) 

Figure  2 shows  a case  with  a fairly 
strong  shock,  the  flow  past  an  NACA  0012 
airfoil  at  Mach  0.750  and  an  angle  of 
attack  of  2*.  In  this  case  the  calculation 
converges  rapidly  after  an  initial  hesita- 
tion. A study  of  plots  of  the  pressure 
distribution  after  each  cycle  chows  that 
the  formation  of  the  shock  is  accompanied 
by  the  appearance  ahead  of  the  shock  of  a 
temporary  overshoot  which  is  subsequently 
suppressed.  The  first  order  accurate 
scheme  obtained  by  setting  c • 1 in 
equation  (11)  was  used  both  in  this  calcu- 
lation and  in  the  calculations  displayed  in 
Figure  1. 

Figure  3 shows  the  initial  convergence 
history  of  the  calculation  of  a flow  con- 
taining two  shocks,  to  illustrate  the 
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spaed  with  which  the  flow  pattern  is 
established.  The  case  is  that  of  a Korn 
airfoil  at  a Mach  number  slightly  below  its 
wC5*?n  c&lcul2tsd  with  the  ssccr.d 

order  accurate  scheme  obtained  by  setting 
t • 1 ui  equation  (11)  . (The  forward  shock 
was  suppressed  when  the  calculation  was 
repeated  with  c » 0).  The  potential  of 
incompressible  flow  generated  by  the 
conformal  mapping  was  used  to  start  the 
calculation,  yielding  the  pressure  distri- 
bution displayed  in  Figure  3(a).  The 
subsequent  plots  show  the  pressure  distri- 
bution after  completing  the  iteration  on 
the  fine  grid  and  before  entering  the 
multigrid  loop  of  each  cycle.  At  the 
beginning  of  the  10th  cycle  the  calcula- 
tion is  essentially  complete.  The  lift 
and  drag  coefficients  CL  » 0.5998  and 
CO  ■ 0.0003  are  identical  to  the  values 
obtained  when  this  calculation  was 
continued  for  50  cycles.  It  appears  that 
it  should  generally  be  possible  to 
calculate  the  flows  likely  to  be 
encountered  in  subsonic  cruising  flight 
with  about  10  cycles  of  the  multigrid 
alternating  direction  scheme. 
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